function fval = int_partial(L, B, a_b, a_s, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, DELTA, CHI_q, CHI_k, CHI_PI, CHI_PI_k, tau_k, Z_K)

sig = S_PARAM;
%mu  = - 0.5 * sig^2;
mu  = M_PARAM;
f   = @(x, a) 1 ./ (x * sig * sqrt(2 * 3.14159265359)) .* exp(- (log(x) - mu - a).^2 / 2 / sig^2)...
   .* (x >= EPS_L_BOUND) .* (x <= EPS_H_BOUND);
int_method = 'auto';
Rel_Tol    = 1e-5;
Abs_Tol    = 1e-6;

fun.l_bound      = @(a) EPS_L_BOUND; % lower bound
fun.h_bound      = @(a) EPS_H_BOUND; % upper bound
fun.Delta        = @(y, x, B, CHI_q, CHI_PI) (1 - tau_k) .* B .* ((1 - THETA - CHI_PI) .* x + THETA .* y) + (1 - DELTA) .* (1 - CHI_q);
fun.threshold    = @(x, L, CHI_PI, CHI_PI_k) L ./ (1 - THETA - CHI_PI - CHI_PI_k) - THETA ./ (1 - THETA - CHI_PI - CHI_PI_k) .* x; 

fun.reall_P      = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K)...
    (fun.Delta(y, x, B, CHI_q, CHI_PI) + CHI_q .* (1 - DELTA) + CHI_PI .* (1 - tau_k) .* B .* x) .* f(y, a_sell) .* f(x, a_buy)...
    .* (Z_K + (1 - DELTA) .* CHI_k + CHI_PI_k .* (1 - tau_k) .* B .* x) ./ fun.Delta(y, x, B, CHI_q, CHI_PI);

int_l = min(max(L / THETA - (1 - THETA - CHI_PI - CHI_PI_k) / THETA *  fun.h_bound(a_b), fun.l_bound(a_s)), fun.h_bound(a_s)); %need an explanation
int_h = min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)); %need an explanation


fval  = integral2(@(eps_t,eps)fun.reall_P(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
            int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
            + integral2(@(eps_t,eps)fun.reall_P(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
            int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b), 'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);
end